#read in DSI from web and convert from HDF to ASCII
#before read in to Stata in read_lights_drought.do



#start cloek
ptm <- proc.time()

library(raster)



for (i in 2000:2011) {
  
  
  infile <- paste("http://files.ntsg.umt.edu/data/NTSG_Products/DSI/Annual/DSI_05deg_", i, ".hdf", sep="")
  outfile <- paste("./data/grid_ascii/drought_", i, ".asc", sep="")

  r <- raster(infile)
  
  r[r > 1e+35] <- NA
  #leave off top 5 degrees to match lights 
  rc <- crop(r, extent(r, 101, 2800, 1, 7200))
  
  writeRaster(rc, filename=outfile, format = "ascii", overwrite=TRUE, NAflag=-9999)

}

# Stop the clock
proc.time() - ptm
